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1  Introduction 

An  automated  technique  has  been  developed  and  evaluated  to  reconstruct  3-D 
binary  images  of  breast  calcifications.  The  reconstruction  algorithm  consists  of 
segmentation,  motion  correction,  correlation  between  views,  3-D  binary  limited- 
view  reconstruction  of  each  calcification,  and  3-D  rendering.  Our  previous  work 
relied  upon  significant  human  intervention  and  judgment  in  producing  the  final 
3-D  image.  In  this  grant,  we  sought  methods  to  automate  these  tasks.  Required 
were  robust  methods  of  identifying,  segmenting  and  correlating  (or  pairing) 
calcifications  between  views. 

The  tasks  of  identifying  and  segmenting  calcifications  have  been  attempted 
on  numerous  occasions.  These  previous  attempts  have  been  used  almost  exclu¬ 
sively  in  computer  aided  diagnosis  (CAD)  systems.  In  such  systems,  the  desire 
is  to  capture  a  sufficient  number  of  calcifications  to  identify  clusters  of  suspicious 
calcifications  for  evaluation  by  a  human  observer.  Significant  effort  is  expended 
on  eliminating  false  positives.  By  imaging  the  breast  at  3  separate  angles  with 
known  spatial  alignment,  we  have  the  advantage  that  true  calcifications  are 
present  in  each  of  the  images,  while  most  spurious  or  non-calcified  signals  are 
only  found  in  one  image.  This  additional  constraint  allows  us  to  segment  more 
calcifications  in  each  image,  admittedly  with  a  higher  fake  positive  rate.  The 
task  of  correlation  between  the  images  quite  naturally  reduces  the  false  positive 
rate  in  the  3-D  image. 

The  work  of  this  grant  is  reviewed  in  this  final  report. 
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2  Body 

2.1  Summary  of  Work  Items 

It  is  useful  to  restate  the  work  items  listed  in  the  original  grant-  They  are  as 
follows: 

Task  1:  Compile  database  of  50  selected  cases  (Months  1-2) 

Task  2:  Manually  identify  and  pair  calcifications  in  database  images  (Months  3-4) 

Task  3:  Evaluate  methods  for  automated  identifications,  segmentation  and  corre¬ 
lation  of  calcifications  (Months  1-24). 

Task  4:  Apply  reconstruction  technique  to  non-calcified  structures  (Months  25-36) 

2.2  Segmentation 

2.2.1  Image  Processing 

This  section  describes  our  algorithm  for  simultaneously  identifying  and  segment¬ 
ing  calcifications  in  individual  images. 

Each  image  is  represented  as  a  1024  x  1024  array  of  12-bit  data,  which 
represents  a  5  cm  by  5  cm  region  in  an  x-ray  projection.  Figure  1  shows  an 
example  of  such  an  image.  While  the  calcifications,  which  can  be  seen  in  the 
upper  left  hand  corner,  are  darker  than  their  immediate  surround,  the  difference 
is  not  so  large  as  the  variation  across  the  image  due  to  the  structures  in  the 
breast  of  much  larger  scale. 

To  remove  these  larger  features  and  thereby  reduce  the  dynamic  range,  a 
median  filtered  image  is  subtracted  from  the  original  image.  That  is,  from 
each  image  pixel  value  is  subtracted  the  median  of  a  31  x  31  region  containing 
the  pixel.  Figure  2  shows  the  image  after  subtraction  of  the  median  filtered 
image  and  after  denoising  the  image.  The  denoising  consists  of  replacing  pixel 
values  which  differ  by  more  than  three  (3)  standard  deviations  from  the  mean, 
calculated  for  the  surrounding  5x5  region.  The  removed  values  are  replaced 
with  the  average  of  the  neighboring  pixels. 

To  further  emphasize  the  high  frequency  components  of  the  image,  a  Lapla- 
cian  operator  is  applied  to  the  image.  At  each  pixel,  a  quadratic  fit  is  made 
to  a  7  x  7  neighborhood,  and  the  coefficients  of  the  un-mixed  quadratic  terms 
are  used  to  estimate  the  Laplacian.  The  result  of  this  operation  is  shown  in 
Figure  3. 
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tions. 


Figure  3:  Digital  image  after  applying  (negative)  Laplacian  operator. 


2,2.2  Segmentation  of  Candidate  Calcifications 

Figure  4  shows  the  Euler  characteristic  for  the  Laplacian  image  as  a  function 
of  threshold.  The  intention  is  that  this  gives  one  a  reasonable  estimate  of  the 
threshold  at  which  to  segment  the  Laplacian  image.  A  discussion  of  the  Euler 
characteristic  and  several  tests  on  simulated  data  is  given  in  section  2.4. 

Figure  5  shows  a  segmentation  obtained  by  including  all  pixels  in  the  Lapla¬ 
cian  image  (Figure  3)  whose  values  are  less  than  the  threshold  at  which  the 
Euler  characteristic  achieves  one  quarter  of  its  maximum  value.  The  cluster  of 
calcifications  in  the  upper  left  hand  corner  is  clearly  visible.  While  a  significant 
number  of  other  pixels  are  also  identified,  even  in  the  segmented  image  most  of 
these  appear  as  a  “dust55  which  a  human  observer  would  understand  is  not  to  be 
consider  calcified  material.  The  intention  of  the  remainder  of  the  algorithm  is  to 
computationally  reject  these  false  positives,  based  upon  the  segmentation  and 
information  in  the  original  image.  First,  however,  note  that  the  segmentation 
is  stable  in  the  sense  that,  if  instead  of  choosing  the  point  at  which  the  Euler 
characteristic  achieves  25%  of  its  maximum,  one  choose  50%,  the  resulting  im¬ 
age  is  quite  similar,  as  can  be  seen  in  Figure  6.  In  order  to  guard  against  any 
significant  variation  in  the  statistics  of  an  image  across  the  image  (for  example, 
in  some  of  our  images  part  of  the  detector  is  directly  exposed  to  the  x-ray  beam) 
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Eu&r  Characteristic  for  Exaision  sat  as  a  function  of  Threshob 


Threshold 

Figure  4;  Euler  characteristic  of  Laplacian  image  as  a  function  of  threshold. 


the  image  is  divided  into  overlapping  250  x  250  pixel  regions  and  the  collection 
of  all  pixels  so  selected  is  used  for  the  segmentation  (the  overlap  corresponds  to 
advancing  the  250  x  250  in  steps  of  one  third  the  size  of  the  window). 

Calcification  candidates  are  then  constructed  as  collections  of  contiguous 
pixels,  using  four-way  continuity.  For  the  x-ray  image  from  which  images  1-6 
have  been  derived,  the  algorithm  (using  a  threshold  of  25%  of  that  at  which  the 
Euler  characteristic  obtains  its  maximum)  gives  1677  candidates. 

2.2,3  Feature  Analysis 

For  each  calcification  candidate,  a  variety  of  features  are  calculated.  The  first 
set  of  features  uses  a  rectangular  region  from  the  denoised  image  containing  the 
pixels  in  the  calcification  candidate  and  extending  a  further  7  pixels  in  each 
direction.  Within  this  region,  all  pixels  that  would  be  segmented  by  a  cut¬ 
off  corresponding  to  35%  of  the  Euler  characteristic  maximum  (this  excludes 
slightly  more  than  the  pixels  in  the  actual  calcification  candidates)  are  used  for 
a  quadratic  fit  which  serves  as  an  estimate  of  the  background  of  the  calcification. 
Thus,  the  pixels  are  identified  by  the  using  the  heuristic  involving  the  Laplacian 
image  and  the  Euler  characteristic  as  a  function  of  threshold,  but  the  region 
being  used  here  is  the  corresponding  region  in  the  de-trended  and  de-noised 
image,  e.g.  Figure  2.  From  this  fit  the  following  features  are  estimated: 

max  depth  The  maximum  difference  between  the  fitted  quadratic  background  and  the 
actual  image  value  at  any  pixel  associated  with  the  calcification  candidate. 


10 


Figure  5:  Segmentation  based  on  Laplacian  image  at  the  point  where  the  Euler 
characteristic  reaches  25%  of  its  maximum* 
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Figure  6:  Segmentation  based  on  Laplacian  image  at  the  point  where  the  Euler 
characteristic  reaches  50%  of  its  maximum. 
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ave  depth  The  average  difference  between  the  fitted  quadratic  background  and  the 
actual  pixel  values. 

Xf  The  value  of  x2  per  degree  of  freedom,  as  a  measure  of  how  well  the 
quadratic  fit  the  background. 

Xcatc  The  sum  of  the  squares  of  the  differences  between  the  values  of  the  pixels 
in  the  candidate  and  the  estimates  of  those  values  from  the  fit  to  the 
background. 

The  distributions  of  maximum  and  average  depth  features  are  shown  in 
figure  7.  The  distributions  of  and  X2caic  3X6  shown  in  figure  8. 

Another  set  of  features  are  calculated  in  an  attempt  to  exclude  false  calcifi¬ 
cation  candidates  which  axe  associated  with  linear  anatomical  features.  Figure  9 
shows  a  linear  feature  in  which  three  presumably  incorrect  candidates  have  been 
identified.  In  this  image,  the  linear  feature  is  roughly  in  line  with  the  long  axes 
of  the  calcifications.  For  each  calcification  candidate,  the  pair  of  pixels  that  are 
separated  by  the  greatest  Euclidean  distance  is  identified,  and  the  line  through 
this  pair  is  considered  to  be  the  long  axis  of  the  calcification.  The  width  of  the 
calcification  is  then  identified  as  twice  the  farthest  distance  of  any  candidate 
pixel  from  the  long  axis.  A  region  7  times  the  length  of  the  calcification  and  8 
times  the  width  Is  then  divided  into  three  regions  parallel  to  the  long  axis  of  the 
calcification,  and  the  average  of  each  region  is  used  as  a  statistic.  The  resulting 
statistics  of  this  procedure  are 

aspect  ratio  The  ratio  of  the  width  of  the  candidate  to  the  length. 


/^online  ^online  5 
Z^+side^+sIde? 

side side  For  these  quantities,  three  rectangular  regions  whose  long  axes  are 
parallel  to  the  long  axis  of  the  candidate  are  used.  The  “on  line55  region 
contains  the  calcification,  while  the  other  regions  He  to  either  side  (the 
sign  is  an  arbitrary  label).  For  each  region,  the  average  pixel  value  fi  and 
the  standard  deviation  a  In  the  set  of  pixel-values  is  calculated. 

The  distributions  of  these  values  are  shown  for  all  candidates  and  selected 
candidates  in  Figures  10  and  11. 

In  an  additional  search  for  large  scale  structures  which  can  produce  false 
positives,  an  additional  set  of  statistics  are  calculated  as  follows.  First,  the 
average  and  population  variance  of  the  pixels  in  the  candidate  calcification  are 
determined.  Second,  the  average  and  population  variance  of  the  pixels  in  an 
approximately  300  x  300  pixel  region  around  the  calcification  candidate  are 
calculated.  Finally,  starting  with  the  calcification,  one  determines  the  largest 
connected  region  such  that  all  of  the  pixels  in  that  region  are  less  than  the  av¬ 
erage  over  the  300  x  300  region  by  one  population  variance.  This  new  region 
is  intended  to  represent  any  large  scale  structure  in  which  the  candidates  are 
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Figure  10:  Aspect  ratio  of  calcification  candidates. 
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Statistics  Related  to  Linear  Features 


Figure  11:  On-line  (^online)  and  off-line  (m+s ide  and  fi-siae  combined).  The 
on-line  statistics  have  been  scaled  by  a  factor  of  2. 


located.  Figure  9  shows  an  example  of  the  region  determined  for  the  three  candi¬ 
dates  in  one  linear  feature.  Regions  similarly  associated  with  other  calcification 
candidates  are  shown  in  Figure  12. 

From  this  one  obtains  the  following  features: 

instruct?  ^struct?  ^struct  The  number,  mean,  and  population  standard  deviation  of  pixels  identified 
as  being  in  a  possible  large  scale  structure  containing  the  calcification 
candidate. 

/xcalc,  <TCaic  The  mean  and  population  variance  of  the  pixels  in  the  calcification  candi¬ 
date. 

preg,  <jreg  The  mean  and  population  variance  of  the  pixels  in  the  300  x  300  region. 

The  distributions  of  Nstruct,  /Struct?  Pc aic?  Mreg  are  shown  in  Figure  13. 

To  remove  false  positives,  calcifications  candidates  are  rejected  for  failing 
any  of  the  following  tests: 

1.  If  the  aspect  ratio  <  0,25,  then  candidate  is  assumed  to  be  a  linear 
anatomic  feature  and  not  a  calcification. 

2.  If  either  of  the  side  averages  is  significantly  greater  than  the  on  line  aver¬ 
age,  then  the  candidate  seems  to  lie  inside  a  linear  feature  and  is  rejected 
as  a  false  positive.  Specifically,  if 


P— side  ""  Monline 
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Figure  12:  Segmentation  of  regions  of  large  scale  structure  which  might  indicate 
that  some  calcification  candidates  are  false  positives. 
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Size  of  PesslW*  Larger  Structure 


Figure  13:  Values  for  the  size  of  possible  larger  structures  containing  calcifi¬ 
cations  (AWuct),  the  average  value  in  such  regions  ^struct,  the  average  digital 
value  in  calcification  candidate  fx caic,  and  the  average  value  in  300  x  300  pixel 
regions  containing  the  candidates  (/ireg).  These  values  are  calculated  from  the 
image  after  unsharp  masking,  so  the  baseline  value  has  been  adjusted  to  2000. 


Figure  14;  Example  of  a  region  in  which  calcification  have  been  identified  and 
segmented* 


or 

M+ side  -  A-line  >  5^a2+side  + a2oriVme 
then  the  candidate  is  rejected* 

3.  If  xlaic  <  25x/?  then  the  candidate  is  rejected.  Note  that  for  uncorrelated 
noise,  this  would  correspond  to  a  cutoff  in  the  signal-to-noise  ratio  of  5. 

4.  If  fiTeg  —  /Zcaic  <  2crreg?  then  the  candidate  is  considered  insufficiently  dense 
to  be  a  true  calcification. 

5.  If  mm  depth  <  &%f  or  max  depth  <  £crreg5  the  candidate  is  similarly 
rejected. 

6.  If  Instruct  >  1000,  the  candidate  is  assumed  to  be  a  false  positive  resulting 
from  a  large  scale  anatomical  structure. 

7.  If  NstTVLCt  is  greater  than  twice  the  number  of  pixels  in  the  candidate  and 

instruct  locale  ^  o* struct  ® calc  then  the  candidate  is  rejected. 

Figure  14  show  a  region  containing  calcifications  that  have  been  segmented. 
Some  examples  of  candidates  which  are  probably  false  positives  are  shown  in 
figure  15. 

2.2,4  Performance 

Having  no  true  “gold  standard*’  available,  we  compared  the  set  of  calcifications 
identified  by  our  algorithm  with  those  identified  by  a  human  observer.  The 


i 


Comparison  of  Human  and  Algorithmic  Observer 
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Figure  16:  Comparison  of  the  number  of  candidates  detected  by  the  algorithm 
with  the  number  identified  by  a  human  observer. 


human  observer  was  a  medical  physicist  familiar  with  mammography.  A  com¬ 
parison  of  the  number  of  calcifications  detected  in  each  image  by  the  human 
observer  with  the  number  algorithmically  detected  is  shown  in  figure  16,  with 
a  correlations  coefficient  of  0.5.  Figure  17  shows  the  fraction  of  calcification  de¬ 
tected  by  the  human  observer  which  were  also  detected  algorithmically,  and  also 
the  fraction  of  the  calcifications  detected  algorithmically  which  were  detected 
by  the  human  observer.  Approximately  one  half  of  the  manually  identified  cal¬ 
cifications  were  identified  algorithmically. 

Retrospectively,  among  the  calcifications  identified  algorithmically  but  not 
identified  by  the  human  observer,  there  was  a  continuum  ranging  from  shadows 
which  were  almost  certaintly  not  calcifications,  through  various  degrees  of  am¬ 
biguity,  to  objects  which  were  almost  certaintly  calcifications  which  had  been 
overlooked  by  the  human  observer.  A  subset  of  the  images  had  been  similarly 
analyzed  by  a  second  human  observer,  also  a  physicist.  Agreement  between 
physicists  was  about  as  good  as  agreement  between  the  algorithm  and  each 
physicist.  We  have  discussed  with  radiologists  the  fact  that  there  seems  to  be 
an  irreducible  degree  of  ambiguity  in  deciding  which  shadows  on  a  mammogram 
represent  calcifications  and  which  ones  don’t.  This  makes  it  difficult  to  “tune” 
an  algorithm  in  a  meaningful  manner. 

2.3  Correlation 

Each  calcification  lies  on  the  geometric  line  from  the  position  of  the  x-ray  focus 
at  the  time  of  acquisition  to  the  position  of  the  shadow  of  the  calcification  in 
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Agreement  Between  Human  and  Algorithm 


Figure  17:  Fraction  of  candidates  identified  by  the  human  or  algorithmic  ob¬ 
server  upon  which  the  two  observers  agree  upon  the  presence  of  a  calcification. 


the  image.  Thus,  for  corresponding  calcifications  in  multiple  images,  the  lines- 
of-sight  connecting  the  x-ray  foci  and  the  projections  must  intersect  at  the 
position  of  the  calcification  in  space.  This  requirement  allows  one  to  identify 
shadows  corresponding  to  the  same  projection  in  multiple  images. 

In  practice,  we  have  found  several  issues  that  confound  this  approach.  One 
issue  is  that  candidates  visible  in  one  image  occasionally  are  not  visible  in  the 
second  image  or  overlap  other  calcifications,  precluding  a  unique  matching  of 
calcifications.  A  second  issue  we  have  found  is  that  projected  shadows  are  often 
not  consistent  with  any  matching  between  images.  This  appears  to  be  due 
to  a  combination  of  patient  motion  and  some  uncertainty  or  instability  in  the 
acquisition  geometry. 

As  a  result,  both  the  matching  problem  and  the  geometry  correction  prob¬ 
lem  must  be  dealt  with  simultaneously.  The  most  general  possible  geometric 
correction  would  introduce  a  large  number  of  parameters  to  fit.  We  have  settled 
on  providing  an  effective  position  adjustment  of  the  data,  consisting  of  a  trans¬ 
lation  of  the  second  off-axis  image  in  the  horizontal  direction  and  translations 
of  the  third  (on-axis)  image  in  both  the  horizontal  and  vertical  directions.  This 
combined  algorithm  is  discussed  in  the  next  subsection. 

2.3.1  Algorithm 

For  a  point  P  and  a  line  1,  let  d(P,  1)  be  the  perpendicular  distance  from  the  point 
to  the  line.  For  a  given  set  of  sight-lines  1%,  I2,  I3,  each  joining  the  position  of 
the  x-ray  focus  to  the  corresponding  shadow  of  a  calcification,  the  calcification 
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should  lie  at  the  common  intersection.  Due  to  experimental  uncertainties,  it 
is  then  reasonable  to  use  as  the  position  of  the  calcification  the  point  which 
minimizes  the  distances  to  the  sight-lines  in  a  least  squares  sense,  ie.  Pcaic  is 
the  point  P  which  minimizes: 

drms{P,h,h,h)  =  vW-P.M  +  CP{P,l2)  +  <P(P,h))  /3.  (1) 

The  minimum  RMS  distance  then  gives  an  measures  of  how  well  the  three  sight¬ 
lines  fit  the  hypothesis  that  they  intersect  at  a  point. 

For  a  given  adjustment  of  the  acquisition  geometry,  it  was  necessary  to  ex¬ 
amine  all  possible  matches  between  shadows  in  each  image.  To  remove  matches 
that  were  clearly  inconsistent,  and  thereby  reduce  the  combinatorial  complexity, 
the  following  preliminary  test  was  performed.  For  each  pair  of  projections  in 
the  off-axis  images,  the  position  of  the  hypothetical  calcification  corresponding 
to  those  two  projections,  Poview,  was  determined  by  minimizing  the  RMS  error 
using  only  the  two  view,  ie.: 

drms(P,  h,h)  =  VW(P,  h)  +  <P(P,  h))  /2.  (2) 

The  point  P2view  was  then  projected  into  the  third  view  and  only  shadows 
within  0.72  cm  of  the  projection  of  the  hypothetical  calcification  where  further 
considered.  For  all  remaining  triples,  the  values  of  dTms  using  the  three  sight¬ 
lines  was  then  computed. 

All  hypothetical  matches  of  shadows  between  the  three  views  were  then 
sorted  according  to  drms.  Any  set  of  three  shadows  in  which  one  of  the  shadows 
had  occurred  in  a  set  of  three  with  a  smaller  dTms  was  then  removed  from  the 
list,  so  that  each  shadow  belonged  at  most  to  one  set  of  three  shadows.  As  a 
measure  of  the  overall  success  in  fitting  all  of  the  calcifications,  we  calculated: 

G°F  =  (ix£,)  '  +  dLs(hM,h)/dtutoS  (3) 

where  dcutoff  =  0.05mm  and  the  sum  runs  over  all  matches  of  sight-lines  which 
passed  the  preliminary  descriminator  as  discussed  in  the  previous  paragraph 
and  drms  dCVL%0^ . 

The  GOF  (Goodness  of  Fit)  score  was  then  used  to  adjust  the  geometric 
acquisition  parameters  to  approximately  compensate  for  .  the  uncontrolled  geo¬ 
metric  effects  discussed  above.  The  geometry  was  varied  by  shifting  the  second 
off-axis  view  up  to  1  cm  in  the  horizontal  direction  from  its  nominal  position, 
and  adjusting  the  third  (on  axis  view)  up  to  1  cm  in  either  the  horizontal  or 
vertical  direction.  Adjustments  were  made  in  steps  of  48  ^m,  corresponding  to 
the  size  of  the  pixels.  The  adjustment  which  gave  the  best  GOF  was  then  taken 
as  the  best  approximate  compensation  for  the  combination  of  patient  movement 
and  geometric  uncertainty. 

2.3.2  Results 

Figure  18  shows  three  views  acquired  for  one  of  our  cases.  Four  of  the  calcifica¬ 
tion  matches  are  shown  by  corresponding  labels  in  the  three  images.  The  first 
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three  matches  (a,  b,  and  c)  are  intuitively  reasonable.  The  fourth  match  (d)  is 
probably  incorrect  in  the  middle  figure,  due  to  the  failure  of  the  segmentation 
algorithm  to  locate  the  corresponding  shadow  in  this  region.  Figure  19  shows 
the  fraction  of  calcification  candidates  automatically  identified  in  each  image 
which  were  matched  with  candidates  in  the  other  two  image.  Figure  20  shows 
the  distribution  of  drms  for  all  calcifications  in  all  images. 

2.4  Euler  Characteristic 

This  section  represents  some  notes  on  the  use  of  the  Euler  characteristic. 

2.4.1  Euler  Characteristic 

For  a  selected  set  of  picture  elements,  the  Euler  characteristic  represents  the 
number  of  contiguous  regions  into  which  that  set  of  elements  can  be  divided 
minus  the  number  of  holes  in  those  objects.  In  particular,  the  Euler  character¬ 
istic  can  be  computed  for  the  excursion  set  (the  set  of  pixels  passing  a  given 
thresholding  condition)  as  a  function  of  threshold.  As  the  Euler  characteristic 
can  be  calculated  in  terms  of  local  data,  the  Euler  characteristic  of  excursion 
sets  as  a  function  of  threshold  can  be  calculated  in  a  single  pass  through  the  im¬ 
age,  making  this  calculation  appealing  for  computational  reasons.  For  stringent 
thresholding  conditions  (so  that  the  excursion  set  is  small)  there  are  generally 
few  holes  in  the  resulting  excursion  set,  so  the  Euler  characteristic  approximates 
the  number  of  objects  that  would  be  in  the  binary  image  after  thresholding.  We 
looked  at  this  with  the  intent  that  by  studying  the  Euler  characteristic  as  a 
function  of  threshold,  one  could  then  choose  a  threshold  which  reveals  mean¬ 
ingful  objects  in  the  resulting  segmentation.  This  is  based  loosely  on  the  fact 
that  human  observers,  when  attempting  to  set  such  a  threshold  manually,  have 
some  success  by  looking  for  a  threshold  at  the  point  at  which  distinct  objects 
seem  to  appear. 

2.4.2  Calculation  of  the  Euler  Characteristic 

The  Euler  characteristic  is  equal  to  the  number  of  regions  minus  the  number 
of  holes.  The  advantage  of  the  Euler  characteristic  is  that  it  can  be  computed 
from  local  data  as 


X  =  #Faces  -  #Edges  +  #Vertices.  *  (4) 

Given  a  set  of  pixels  S  in  a  square  grid,  two  methods  suggest  themselves  for 
computing  the  Euler  characteristic.  For  the  first  method,  Xa(S),  one  considers 
each  pixel  to  represent  a  square  region  of  the  plane.  If  the  pixel  is  in  the  set  S, 
then  that  pixel  is  included  among  the  set  of  Faces,  its  four  edges  are  included 
among  the  set  of  Edges,  and  its  four  vertices  axe  counted  among  the  set  of 
Vertices.  As  each  edge  is  shared  by  two  pixels  (except  at  boundaries)  and  each 
vertex  is  shared  by  four  pixels  (except  at  boundaries)  one  must  take  care  to 
not  count  edges  or  vertices  multiple  times.  A  second  method,  Xb(S ),  results 
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Figure  18:  Segmentation  in  3  views.  Four  three-way  matches  are  shown,  with 
root-mean-square  distances  to  lines  of  (a)  0.002  cm,  (b)  0.002  cm  (c)  0.003  cm 
and  (d)  0.03  cm. 
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Figure  19:  Percent  of  selected  calcification  candidates  automatically  matched. 
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Figure  20:  Root-mean-square  distance  to  lines-of-sight  from  fitted  positions 
calcifications 
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Figure  21:  This  set  of  pixels  can  count  as  one  or  two  objects:  xa  —  1  and 
Xb  =  2. 


from  considering  each  pixel  value  to  be  a  sample  point  at  the  center  of  the 
pixel.  The  number  of  vertices  is  then  given  by  the  number  of  pixels  in  the  set  S, 
the  number  of  edges  is  given  by  the  set  of  horizontal  or  vertical  edges  between 
adjacent  pixels  in  S,  and  the  number  of  Faces  is  given  by  sets  of  four  pixels  with 
coordinates  (»,  j),  ( i,j  + 1),  (i  +  1  ,j)  and  (i  + 1  ,j  +  1)  all  of  which  are  in  5. 

The  values  of  xa  and  xb  are  not  identical,  as  seen  in  Figure  21.  The  set 
of  pixels  S  in  this  image  could  be  reasonably  grouped  into  one  or  two  objects. 
According  to  the  counting  rules  for  xa,  there  are  a  total  of  two  faces,  eight 
edges,  and  seven  vertices  (noting  not  to  double  count  the  shared  corner),  giving 
Xa(S)  =  1.  According  to  the  xb  rules,  there  are  no  faces  or  edges,  but  there  are 
two  vertices  (corresponding  to  the  sample  points  at  the  centers  of  the  pixels)  so 
that  xb(S)  =  2. 

The  counting  rules  are  dual  in  the  sense  that  if  the  entire  square  grid  is 
divided  into  two  disjoint  sets  S  and  S',  then 

Xa{S)  +  xb(S')  «  0,  (5) 

i.e.  every  object  in  S  is  a  hole  in  S'  and  vice-versa.  The  equality  can  not  be 
exact  due  to  the  presence  of  a  boundary,  as  connected  regions  of  one  set  which 
reach  the  boundary  axe  not  counted  as  holes  in  the  other  set.  In  order  to  remove 
the  issue  of  boundary  conditions,  for  the  remainder  of  this  paragraph  we  will 
assume  periodic  boundary  conditions.  The  equality  in  equation  5  is  exact  under 
these  boundary  conditions,  so  that  the  approximation  in  equation  5  under  the 
boundary  actually  implemented  is  simply  related  to  ignoring  boundary  effects. 
Now  to  understand  equation  5,  observe  that  every  pixel  is  either  in  S  or  S'.  If 
it  is  in  S,  then  in  contributes  a  face  according  to  the  xa  counting  rules,  while 
if  it  is  in  S'  then  it  contributes  a  vertex  according  to  the  xb  counting  rules,  so 
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that 


#Faces  A(S)  +  #Vertices  B(ff)  =  N  (6) 

where  N  is  the  total  number  of  pixels  in  the  grid.  Further,  given  two  adjacent 
pixels  (either  horizontal  or  vertical)  then,  it  is  either  the  case  that  at  least  one 
of  the  pixels  is  in  S,  so  that  the  common  edge  between  the  pixels  contributes 
to  the  number  of  edges  according  to  the  xa  rules,  or  both  pixels  axe  in  S',  in 
which  case  the  edge  between  the  sample  points  contributes  to  the  number  of 
edges  according  to  the  xb  rules,  and  as  these  cases  are  mutually  exclusive  one 
has 

#Edg esA(S)  +  #Edges  B(S>)  =  JVedge  (7) 

where  Ne dge  is  the  total  number  of  edges.  Similarly,  every  corner  of  every  pixel 
is  either  the  corner  of  a  pixel  in  S  or  is  the  corner  of  four  pixels  with  sample 
points  in  S',  and  again  as  these  cases  are  mutually  exclusive  one  has 

#  Vertices^  (S)  +  #FacesB(S')  =  NveTtices  (8) 

where  Arvertices  is  the  number  of  corners  of  all  pixels  in  the  image.  Combining 
equations  6,  7,  and  8  and  noting  that,  with  periodic  boundary  conditions  N  — 
iV€dge  +  ^vertices  =  0,  one  obtains  the  result  given  by  equation  5. 

The  ambiguity  in  choosing  between  the  counting  rules  corresponding  to  xa 
and  xb  produces  several  odd  features.  Note,  for  example,  that  if  the  set  of  pixels 
shown  in  gray  in  Figure  21  axe  take  as  the  set  of  pixels  of  interest,  then  one  object 
is  visible,  but  if  the  two  pixels  axe  taken  as  the  excluded  set,  then  there  are  two 
holes  in  the  region.  This  asymmetry  produces  some  counter-intuitive  results. 
For  example,  if  the  random  process  generating  the  field  of  pixels  values  produces 
values  symmetrically  relative  to  some  mean  value  (as,  for  example,  a  Gaussian 
random  field)  then  one  would  intuitively  expect  that  the  Euler  characteristic  for 
the  excursion  set  as  a  function  of  threshold  would  be  symmetric  with  respect 
to  the  mean  value  of  the  field.  However,  because  of  situations  like  the  above, 
configurations  like  that  in  Figure  21  carry  more  weight  as  “holes”  than  as  pixels 
that  are  part  of  the  segmentation,  so  that  the  curve  is  not  symmetric.  Instead, 
the  Euler  characteristic  tends  to  be  negative  at  a  threshold  corresponding  to 
the  image  mean  and  the  magnitude  of  the  most  negative  value  is  greater  than 
the  magnitude  of  the  highest  value. 

It  is  interesting  to  note  that  the  following  change  of  rules  removes  this  ambi¬ 
guity.  Consider  the  pixel  values  localized  on  “sample  points55 ,  each  sample  point 
being  a  possible  vertex  in  the  selected  set.  As  before,  let  edges  be  lines  between 
adjacent  vertices.  However,  choose  the  lines  between  adjacent  vertices  in  such 
a  way  that  the  resulting  set  of  vertices  and  line  segments  results  in  dividing  the 
plane  up  into  triangles  (requiring,  of  course,  that  two  distinct  segments  meet 
only  at  vertices  corresponding  to  sample  points  and  that  no  point  is  inside  two 
distinct  triangles).  Now  use  the  resulting  triangles  as  candidates  for  faces,  the 
triangle  serving  as  a  face  if  and  only  if  all  three  vertices  are  in  the  set  5.  Under 
these  rules,  there  are  no  ambiguities  analogous  to  that  in  Figure  21,  and  one 
can  prove  that  (again  requiring  periodic  boundary  conditions) 

x(S)+x(S')  =  0  (9) 
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where  the  Euler  characteristics  for  both  sets  are  computed  by  the  same  rules, 
as  described  earlier  in  the  paragraph. 

2.4.3  Random  Field  Examples 

Figure  22  shows  each  component  of  the  Euler  characteristic  and  the  Euler  char¬ 
acteristic  itself  (according  to  the  xa  rules )  as  a  function  of  threshold,  normalized 
by  the  number  of  pixels  in  the  image.  The  simulation  in  Figure  22  corresponds 
to  white  Gaussian  noise  on  a  2048  x  2048  grid.  As  is  well  known [1],  the  curves 
in  Figure  22  can  be  calculated  from  a  simple  theoretical  argument.  Given  that 

f(*<0  =  l(l  +  ert(i))  (10) 

is  the  probability  of  a  given  pixel  having  value  x  less  than  or  equal  to  a  thresh¬ 
old  t,  one  can  treat  the  inclusion  of  each  possible  face,  edge,  or  vertex  as  a 
independent  random  event  and  obtain 

Faces  (i)/iV  =  P(x  <  t)  (11) 

Edges(t)/JV  =  1  -  (1  -  P{x  <  t)f  (12) 

Vertices (t)/N  =  1  -  (1  -  P(x  <  i))4  (13) 

where  N  is  the  number  of  pixels  in  the  image.  When  so  calculated,  agreement 

was  found  to  be  better  than  1%. 

Figure  23  shows  the  Euler  characteristic  as  a  function  of  threshold  for  four 
different  noise  sources.  In  each  case,  the  threshold  has  been  rescaled  so  that 
the  horizontal  axis  is  in  units  of  pixel— standard— deviations  relative  to  the  image 
mean.  The  Gaussian  noise  is  the  same  as  in  figure  22.  The  negative  exponential 
and  “power-of-cosine”  curves  correspond  to  two  cases  of  noise  which  decreases 
with  frequency,  though  in  the  first  case  the  noise  power  spectrum  is  not  smooth 
at  zero  frequency  and  in  the  second  case  it  is.  The  “CT”  noise  corresponds  to 
noise  which  ramps  up  with  frequency,  somewhat  like  CT  noise.  In  the  “pow- 
cos”case  the  noise  results  in  very  strong  correlations  and  a  few  large  regions,  so 
the  vertical  axis  has  been  rescaled. 

2.4.4  Discussion  and  Summary  of  Scientific  Results 

In  this  grant,  we  generated  a  manually  segmented  and  paired  dataset  of  110 
patients  images,  which  we  have  used  as  a  “gold  standard”  in  the  evaluation  of 
computer  algorithms  for  identifying,  segmenting  and  correlating  calcifications. 
We  have  been  able  to  develop  two  separate  computer  algorithms,  one  for  seg¬ 
mentation  of  the  regions,  the  other  for  correlationg  calcifications  between  the 
images.  Both  are  quite  robust.  There  are  a  number  of  significant  findings  from 
this  work  that  will  be  published.  First,  the  use  of  the  Euler  characteristic  to 
determine  connectivity  in  an  automated  fashion  is  unique.  Secondly,  the  simul¬ 
taneous  correction  of  patient  motion  and  the  determination  of  correspondence 
between  the  views  is  unique. 
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Euler  characterstic  vs.  Threshold 

{for  white  noise,  normalized  to  number  of  pixels) 


Figure  22:  %A  and  components  as  a  function  of  threshold  for  white  noise.  Values 
have  been  divided  by  the  number  of  pixels  in  the  region.  The  horizontal  axis  is 
in  units  of  standard  deviations  of  the  pixel  values. 


3  Key  Research  Accomplishments 

The  following  is  a  list  of  key  research  accomplishments  resulting  from  this  work: 

1*  Developed  a  database  of  110  biopsy-proven  cases,  with  3  digital  images  of 
each  case 

2.  Developed  a  set  of  segmented  images  from  each  of  110  cases.  In  these 
data,  each  calcification  from  all  three  views  of  each  patient  was  manually 
identified,  and  semi-automatically  segmented. 

3.  Developed  a  set  of  manually  determined  correspondences. 

4.  These  datasets  were  used  to  develop  an  automatic  identification  and  seg¬ 
mentation  algorithm  that  tested  each  point  in  an  image  as  a  potential 
seed  point  and  then  tested  each  resultant  segmented  region  for  validity  as 
a  potential  calcification.  A  key  feature  of  this  algorithm  was  the  use  of 
the  Euler  characteristic  to  determine  connectivity. 

5.  The  above  datasets  were  also  used  to  develop  an  automatic  correspon¬ 
dence  algorithm.  The  algorithm  used  a  freighted  summation  that  allowed 
us  to  simultaneously  correct  for  patient  motion  and  determine  optimal 
correspondence. 
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5  Conclusions 

In  conclusion,  we  have  developed  automated  algorithms  for  identifying,  seg¬ 
menting,  and  correlating  calcifications  in  3-D,  using  3  source  images  acquired  at 
15  degree  increments.  The  algorithms  have  been  tested  with  previously  acquired 
clinical  data,  which  was  arranged  into  a  database,  and  was  analyzed  by  human 
observers  for  the  purpose  of  developing  a  gold  standard  for  the  reconstructions. 
The  algorithms  have  worked  very  well.  The  use  of  the  Euler  characteristic  for 
connectivity  analysis  and  the  simultaneous  correction  of  image  correlation  and 
image  motion  are  particularly  noteworthy  accomplishments. 
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